{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 支持向量机（SVM）\n",
    "### 1. 基本概念\n",
    "\n",
    "**支持向量：** 在线性可分的情况下，训练数据集中与分类超平面距离最近的样本点的实例。\n",
    "**间隔**\n",
    "**间隔边界**\n",
    "**函数间隔：** 对于给定的训练数据集$T$和超平面$(w,b)$，定义超平面和样本点$(x_i,y_i)$的函数间隔为：\n",
    "$$\n",
    "\\hat{\\gamma}_i=y_i(w\\cdot{x_i}+b) \\tag{1}\n",
    "$$\n",
    "所有样本点关于超平面的函数间隔可以表示分类预测正确性和确信度：\n",
    "$$\n",
    "\\hat{\\gamma}= \\mathop{min}\\limits_{i=1,...,N}\\hat{\\gamma}_i \\tag{2}\n",
    "$$\n",
    "该间隔当$w$和$b$等比例扩大时，间隔也扩大，所以引入**几何间隔**，使间隔度量确定。\n",
    "\n",
    "**几何间隔：** 超平面$(w,b)$关于样本点$(x_i,y_i)$的几何间隔一般使实例到超平面的带符号距离，当样本点被超平面分类正确时就是实例点到超平面的距离。\n",
    "$$\n",
    "{\\gamma}_i=y_i\\left(\\frac{w}{||w||}\\cdot x_i + \\frac{b}{||w||}\\right) \\tag{3}\n",
    "$$\n",
    "\n",
    "$$\n",
    "{\\gamma}= \\mathop{min}\\limits_{i=1,...,N}{\\gamma}_i \\tag{4}\n",
    "$$\n",
    "\n",
    "由上可以函数间隔和几何间隔的关系：\n",
    "$$\n",
    "\\gamma=\\frac{\\hat{\\gamma}}{||w||}\n",
    "$$\n",
    "**硬间隔最大化：** 对训练数据找到几何间隔最大化的超平面意味着：以充分大的确信度对训练数据进行分类。\n",
    "$$\n",
    "\\begin{aligned}\n",
    "  &\\max\\limits_{w,b} \\gamma \\\\\n",
    "  &s.t.\\  \\ \\  y_i\\left(\\frac{w}{||w||}\\cdot x_i + \\frac{b}{||w||}\\right) \\ge \\gamma    \n",
    "\\end{aligned}\\tag{5}\n",
    "$$\n",
    "由函数间隔和几何间隔的关系可以将上式改为如下凸二次规划问题：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "  &\\mathop{min}\\limits_{w,b} \\quad \\frac{1}{2}||w||^2  \\\\\n",
    "  &s.t. \\quad y_i(w\\cdot x_i+b)-1 \\ge 0, \\quad i = 1,...,N\n",
    "\\end{aligned}\\tag{6} \n",
    "$$\n",
    "\n",
    "**硬间隔最大化对偶形式：** 支持向量机学习的基本思想是求解找到能够正确划分训练数据集并且几何间隔最大的分离超平面。\n",
    "![](img/![](img/2020-03-18-16-23-36.png).png)\n",
    "\n",
    "---\n",
    "\n",
    "### 3. 软间隔最大化 \n",
    "软间隔最大化的凸二次规划形式为：\n",
    "![](img/2020-03-18-16-35-54.png)\n",
    "\n",
    "更改为对偶形式后算法为：\n",
    "![](img/2020-03-18-16-24-34.png)\n",
    "![](img/2020-03-18-16-25-06.png)\n",
    "### 4. 非线性支持向量机\n",
    "![](img/2020-03-18-16-25-50.png)\n",
    "![](img/2020-03-18-16-26-32.png)\n",
    "![](img/2020-03-18-16-27-10.png)\n",
    "### 5. SMO算法\n",
    "SMO算法包括两个部分：求解两个变量的二次规划的解析方法、选择变量的启发式方法。\n",
    "**两变量二次规划求解**\n",
    "![](img/2020-03-18-17-20-29.png)\n",
    "![](img/2020-03-18-17-21-00.png)\n",
    "![](img/2020-03-18-17-23-43.png)\n",
    "<font color='red'>\n",
    "**算法（SMO算法）**\n",
    "\n",
    "&emsp;&emsp; 输入：训练数据集$T= {(x_1,y_1 )，(x_2,y_2 )，…，(x_N,y_N )}$，其中，$x_i\\in \\mathcal{X}= R^n$，$y_i \\in \\mathcal{Y} = \\{-1,+1\\}, i=1,2,…,N$，精度$\\epsilon$:\n",
    "\n",
    "&emsp;&emsp; 输出：近似解$\\hat{\\alpha}$\n",
    "\n",
    "&emsp;&emsp; (1) 取初值$\\alpha^{(0)}=0$，令$k=0$；\n",
    "\n",
    "&emsp;&emsp; (2) 选取优化变量$\\alpha_1^{(k)}, \\alpha_2^{(k)}$，解析求解两个变量的最优化问题(7.101)~(7.103)，求得最优解$\\alpha_1^{(k+1)}, \\alpha_2^{(k+1)}$;\n",
    "\n",
    "&emsp;&emsp;&emsp; a) 外循环选择并更新$\\alpha_1^{(k)}$，内循环选择并更新$\\alpha_2^{(k)}$;\n",
    "\n",
    "&emsp;&emsp;&emsp; b) 更新$\\alpha$为$\\alpha^{(k+1)}$，并且截断（公式7.108）;\n",
    "\n",
    "![](img/2020-03-19-14-36-15.png)\n",
    "![](img/2020-03-19-14-36-47.png)\n",
    "![](img/2020-03-19-14-24-57.png)\n",
    "&emsp;&emsp;&emsp; c) 更新$b$和所有样本的误差$E_i$;\n",
    "\n",
    "![](img/2020-03-19-14-23-32.png)\n",
    "\n",
    "&emsp;&emsp;（3）若在精度$\\epsilon$范围内满足停机条件\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "  &\\sum_{i=1}^{N}\\alpha_i y_i = 0, \\\\\n",
    "  &0\\le\\alpha_i \\le C, i = 1,2...,N \\\\ \n",
    "  &y_i\\cdot g(x_i) \\left\\{ \\begin{aligned}\n",
    "    &\\ge 1, &\\{ x_i|\\alpha_i=0\\} \\\\\n",
    "    &= 1, &\\{ x_i|0<\\alpha_i<C\\} \\\\\n",
    "    &\\le 1, &\\{ x_i|\\alpha_i=C\\}\n",
    "  \\end{aligned}\\right.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "&emsp;&emsp; (3) 否则转到(2)\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.datasets import load_iris\n",
    "from sklearn.model_selection import train_test_split\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 生成数据\n",
    "该数据有两种，一种是使用sklearn自带的鸢尾花数据集，该数据集比较简单，可以快速验证效果\n",
    "\n",
    "另一种是mnist手写字符数据集"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 使用sklearn的鸢尾花数据集，生成数据（前两种label）\n",
    "def gen_data(data_size=100,visualization=False):\n",
    "    iris_data = load_iris()\n",
    "    df = pd.DataFrame(iris_data.data, columns=iris_data.feature_names)\n",
    "    df['label']=iris_data.target\n",
    "    df.columns=['sepal length', 'sepal width', 'petal length', 'petal width', 'label']\n",
    "    data = np.array(df.iloc[:100,[0,1,-1]])\n",
    "    for i in data:\n",
    "        if i[-1] == 0:\n",
    "            i[-1] = int(-1)\n",
    "        elif i[-1] == 1:\n",
    "            i[-1] = int(1)\n",
    "    # print(data)\n",
    "    X = data[:, :2]\n",
    "    y = data[:,-1]\n",
    "    # print(X,y)\n",
    "    if visualization == True:\n",
    "        X1,X2 = X[:50,:],X[50:,:]\n",
    "        label1, label2 = '-1','1'\n",
    "        X_visualization(X1,X2,label1,label2)\n",
    "    return X,y\n",
    "\n",
    "def load_mnist(train_file, test_file):\n",
    "    X_train=[]\n",
    "    X_test=[]\n",
    "    y_train=[]\n",
    "    y_test=[]\n",
    "    with open(train_file,\"r\") as train_f:\n",
    "        data = train_f.readlines()\n",
    "        for line in data:\n",
    "            line = line.strip().split(\",\")\n",
    "            # print(line)\n",
    "            if (int(line[0]) == 0):\n",
    "                y_train.append(1)\n",
    "            else:\n",
    "                y_train.append(-1)\n",
    "            X_train.append([int(i)/255 for i in line[1:]])\n",
    "    with open(test_file,\"r\") as test_f:\n",
    "        data = test_f.readlines()\n",
    "        for line in data:\n",
    "            line = line.strip().split(\",\")\n",
    "            if (int(line[0]) == 0):\n",
    "                y_test.append(1)\n",
    "            else:\n",
    "                y_test.append(-1)\n",
    "            X_test.append([int(i)/255 for i in line[1:]])\n",
    "    # print(len(X_train[0]))\n",
    "    X_train = np.array(X_train)\n",
    "    y_train = np.array(y_train)\n",
    "    X_test = np.array(X_test)\n",
    "    y_test = np.array(y_test)\n",
    "    return X_train,X_test,y_train,y_test\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 可视化数据\n",
    "def X_visualization(X1,X2,label1,label2, show_classification_result=False,w1=0,w2=0,b=0):\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.scatter(X1[:,0],X1[:,1])\n",
    "    ax.scatter(X2[:,0],X2[:,1])\n",
    "    if show_classification_result==True:\n",
    "        x = np.arange(4.0,7.0,0.1)\n",
    "        y = (-b - w1*x) /w2\n",
    "        ax.plot(x,y)\n",
    "        plt.legend(('result',label1, label2),loc='upper right')\n",
    "    else:\n",
    "        plt.legend((label1, label2),loc='upper right')\n",
    "    plt.xlabel('length(cm)')\n",
    "    plt.ylabel('length(cm)')\n",
    "    plt.title(\"iris data\")\n",
    "    plt.show()\n",
    "\n",
    "# 展示分类超平面\n",
    "def show_result(X,w1=1,w2=1,b=1):\n",
    "    X1,X2 = X[:50,:],X[50:,:]\n",
    "    label1, label2 = '-1','1'\n",
    "    X_visualization(X1,X2,label1,label2,show_classification_result=True,w1=w1,w2=w2,b=b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X+cXHV97/HXu0skq/xYA6FINjRUfdBbgYJEkJuWW8VKUYpcVJCrrSkot31ooReLV7w+QCktKhYUsbWgV1QQjUhTVAR/QOov4LohMRED96rFZgNKCCYBDRLC5/5xzk42w+zunJ35zpxz5v18PPaRmTPfPfs5Z2A/O+d8Pt+vIgIzMzOA3+h3AGZmVh5OCmZm1uCkYGZmDU4KZmbW4KRgZmYNTgpmZtbgpGADRdI9kv5wmtffKeljXfg510i6uNP9mPXabv0OwKyXIuIFM7z+972KZYKkFcC1EdFxMjLrlD8pmOUk+Y8kG3hOCjZQJN0v6WX543dLukHStZK2Akvzbdfmr8/NX9skabOk70n6zSn2e4SkuyU9KulzwNxJrz1b0pckbZT0i/zxaP7a3wF/AFwp6TFJV+bbPyRpvaStklZK+oO0Z8Ys46Rgg+5VwA3ACHBd02tvBPYGFgL7AH8BbGvegaRnAMuBTwPzgM8Dr5405DeATwC/BRyY7+NKgIj4X8C3gLdGxB4R8db8e74HHJ7v7zPA5yXNxSwxJwUbdHdExPKIeCoimn/hbydLBs+LiB0RsTIitrbYx4uBOcAHI2J7RNxA9ksdgIjYFBFfiIhfRcSjwN8B/2W6oCLi2vz7noyIfwB2Bw7u4DjN2uKkYINu/TSvfRq4FfispAckvV/SnBbjDgA2xK6zS/504oGkZ0r6Z0k/zS9TfRMYkTQ01Q+W9DZJ6yRtkbSZ7BPLvkUOzGw2nBRs0E05TXD+V/97IuJ3gf8MnAj8WYuhDwILJGnStgMnPX4b2V/5R0fEXsCx+faJ8bvEkN8/+J/AqcCzI2IE2DJpvFkyTgpmU5D0EkmH5n/RbyW7nLSjxdA7gCeBsyXtJukU4KhJr+9Jdh9hs6R5wIVN3/9z4Lebxj8JbAR2k3QBsFc3jslsJk4KZlPbn+wm9FZgHfBvwLXNgyLiCeAUYCnwC+A04MZJQz4IDAMPA3cCtzTt4kPAa/LKpCvILll9Bfi/ZJehHmf6y1xmXSMvsmNmZhP8ScHMzBqcFMzMrMFJwczMGpwUzMysIfkEYHk53xhZc8+JTa8tBS4FNuSbrpxppsh99903Fi1alCBSM7P6Wrly5cMRMX+mcb2YFfIcsnK+qeqsPzdpvpcZLVq0iLGxsa4EZmY2KCT9dOZRiS8f5TNBvhLwPPFmZhWQ+p7CB4G3A09NM+bVktbkUxgvbDVA0lmSxiSNbdy4MUmgZmaWMClIOhF4KCJWTjPsi8CiiDgM+DrwyVaDIuKqiFgcEYvnz5/xkpiZmc1SynsKS4CTJL2CbMGRvSRdGxFvmBgQEZsmjb8aeF/CeMzMCtm+fTvj4+M8/vjj/Q6lbXPnzmV0dJQ5c1pN6DuzZEkhIs4HzgfIF0r/m8kJId/+nIh4MH96EtkNaTOzUhgfH2fPPfdk0aJF7DoJbjlFBJs2bWJ8fJyDDjpoVvvoeZ+CpIsknZQ/PVvSPZK+D5xNNqGYmVkpPP744+yzzz6VSAgAkthnn306+mTTk4XKI2IFsCJ/fMGk7Y1PE2Z1s3zVBi699T4e2LyNA0aGOe/4gzn5iAX9DssKqkpCmNBpvD1JCmaDZvmqDZx/41q2bc+WX9iweRvn37gWwInBSs3TXJglcOmt9zUSwoRt23dw6a339Skiq7p7772XY445ht13350PfOADyX6OPymYJfDA5m2FtpvNZN68eVxxxRUsX7486c/xJwWzBA4YGS603eph+aoNLHnvbRz0ji+z5L23sXzVhpm/qU377bcfL3rRi2ZdatouJwWzBM47/mCG5wztsm14zhDnHX9wnyKy1CbuI23YvI1g532kbiaGXnBSMEvg5CMWcMkph7JgZBgBC0aGueSUQ32Tucbqch/J9xTMEjn5iAVOAgMkxX2kj3zkI1x99dUA3HzzzRxwwAGz3le7/EnBzKwLUtxHestb3sLq1atZvXp1TxICOCmYmXVF6vtIP/vZzxgdHeWyyy7j4osvZnR0lK1bt3Zl35P58pGZWRdMXCpM1cW+//77Mz4+3pV9TcdJwcysS+pwH8mXj8zMrMFJwczMGpwUzMyswUnBzMwanBTMzKzBScEGXspJzMw6dcYZZ7DffvtxyCGH9OTnOSnYQKvLJGZWX0uXLuWWW27p2c9zUrCBVpdJzKwk1iyDyw+Bd49k/65Z1vEujz32WObNm9eF4Nrj5jUbaF4Mx7pmzTL44tmwPf9vZ8v67DnAYaf2L66C/EnBBpoXw7Gu+cZFOxPChO3bsu0V4qRgA82L4VjXbJliXqKptpeULx/ZQEs9iZkNkL1Hs0tGrbZXiJOCDbw6TGJmJXDcBbveUwCYM5xt78Dpp5/OihUrePjhhxkdHeU973kPZ555ZofBTs1Jwfpm+aoN/gvd6mPiZvI3LsouGe09miWEDm8yX3/99V0Irn1OCtYXE/0BE+WgE/0BgBODVddhp1aq0qgV32i2vnB/gFk5OSlYX7g/wKoiIvodQiGdxuukYH3h/gCrgrlz57Jp06bKJIaIYNOmTcydO3fW+/A9BeuL844/eJd7CuD+ACuf0dFRxsfH2bhxY79DadvcuXMZHZ19GayTgvWF+wOsCubMmcNBBx3U7zB6yknB+sb9AWblkzwpSBoCxoANEXFi02u7A58CjgQ2AadFxP2pYzIrG/dsWFn04kbzOcC6KV47E/hFRDwPuBx4Xw/iMSsVr+lgZZI0KUgaBV4JfGyKIa8CPpk/vgE4TpJSxmRWNu7ZsDJJ/Unhg8DbgaemeH0BsB4gIp4EtgD7NA+SdJakMUljVaoCMGuHezasTJIlBUknAg9FxMrphrXY9rSC4Ii4KiIWR8Ti+fPndy1GszJwz4aVScpPCkuAkyTdD3wWeKmka5vGjAMLASTtBuwNPJIwJrPS8ZoOVibJkkJEnB8RoxGxCHgdcFtEvKFp2E3AG/PHr8nHVKN10KxLTj5iAZeccigLRoYRsGBkmEtOOdTVR9YXPe9TkHQRMBYRNwEfBz4t6UdknxBe1+t4zMrAPRtWFj1JChGxAliRP75g0vbHgdf2IgYbHO9avpbr71rPjgiGJE4/eiEXn3xov8MyqwR3NFutvGv5Wq698z8az3dENJ47MZjNzLOkWq1cf1eLNXKn2W5mu3JSsFrZMUWdwlTbzWxXTgpWK0NTNMRPtd3MduWkYLVy+tELC203s135RrPVysTNZFcfmc2OqtYrtnjx4hgbG+t3GGZmlSJpZUQsnmmcPylYV73+6jv4zo93zlSy5LnzuO7Nx/Qxov7xGglWRb6nYF3TnBAAvvPjR3j91Xf0KaL+8RoJVlVOCtY1zQlhpu115jUSrKqcFMwS8BoJVlVOCmYJeI0EqyonBeuaJc+dV2h7nXmNBKsqJwXrmuvefMzTEsCgVh95jQSrKvcpmJkNAPcpWF+kqs0vsl/3B5jNnpOCdc1Ebf5EKeZEbT7Q0S/lIvtNFYPZoPA9BeuaVLX5Rfbr/gCzzjgpWNekqs0vsl/3B5h1xknBuiZVbX6R/bo/wKwzTgrWNalq84vs1/0BZp3xjWbrmokbud2u/Cmy31QxmA0K9ymYmQ0A9ymUVBVr6KsYs5nNjpNCD1Wxhr6KMZvZ7PlGcw9VsYa+ijGb2ew5KfRQFWvoqxizmc2ek0IPVbGGvooxm9nsOSn0UBVr6KsYs5nNnm8091AVa+irGLOZzZ77FMzMBkDf+xQkzQW+Ceye/5wbIuLCpjFLgUuBDfmmKyPiY6listl51/K1XH/XenZEMCRx+tELufjkQzseW5b+h7LEYVYGKS8f/Rp4aUQ8JmkO8G1JX4mIO5vGfS4i3powDuvAu5av5do7/6PxfEdE43nzL/siY8vS/1CWOMzKItmN5sg8lj+dk39V61qVcf1d69veXmRsWfofyhKHWVkUSgqSniVpaOaRjfFDklYDDwFfi4i7Wgx7taQ1km6QtHCK/ZwlaUzS2MaNG4uEbB3aMcU9p1bbi4wtS/9DWeIwK4tpk4Kk35D03yR9WdJDwL3Ag5LukXSppOdP9/0RsSMiDgdGgaMkHdI05IvAoog4DPg68Mkp9nNVRCyOiMXz589v99isC4aktrcXGVuW/oeyxGFWFjN9UrgdeC5wPrB/RCyMiP2APwDuBN4r6Q0z/ZCI2AysAP64afumiPh1/vRq4Mhi4Vtqpx/d8sNby+1Fxpal/6EscZiVxUw3ml8WEdubN0bEI8AXgC/kN5GfRtJ8YHtEbJY0DLwMeF/TmOdExIP505OAdUUPwNKauEHcTkVRkbFl6X8oSxxmZdF2n4KkZwMLmZRIIuLuacYfRnY5aIjsE8myiLhI0kXAWETcJOkSsmTwJPAI8JcRce90cbhPwcysuHb7FNpKCpL+FlgK/JidFUQRES/tJMjZqHpSSFUTX6Q/IOW+ixxfFc9F5axZBt+4CLaMw96jcNwFcNip/Y7K+qDbzWunAs+NiCc6C2uwpaqJL9IfkHLfRY6viueictYsgy+eDdvzSqot67Pn4MRgU2q3JPUHwEjKQAZBqpr4Iv0BKfdd5PiqeC4q5xsX7UwIE7Zvy7abTaHdTwqXAKsk/YCsUxmAiDgpSVQ1laomvkh/QMp9Fzm+Kp6LytkyXmy7Ge0nhU+SVQ6tBZ5KF069HTAyzIYWv/Q6rYkfklr+0puqbyDVvoscXxXPReXsPZpdMmq13WwK7V4+ejgiroiI2yPi3ya+kkZWQ6lq4ov0B6Tcd5Hjq+K5qJzjLoA5TUl2znC23WwK7X5SWJmXj97ErpePpixJtadLVRNfpD8g5b6LHF8Vz0XlTNxMdvWRFdBuSertLTa7JNXMrCK6WpIaES/pPCSrqjL0HljFuV+iMtq6pyDp7yWNTHr+bEkXpwvLymKin2DD5m0EO/sJlq/a0NFYGyAT/RJb1gOxs19izbJ+R2YttHuj+YR8UjsAIuIXwCvShGRlUobeA6s490tUSrtJYUjS7hNP8gnudp9mvNVEGXoPrOLcL1Ep7SaFa4FvSDpT0hnA15hi7QOrlyLrDXhtAmtpqr4I90uUUltJISLeD1wM/CfgBcDf5tus5srQe2AV536JSpm2+kiSIq9ZjYhbgFumG2P1U4beA6s490tUyrR9CpJWkC2m868R8R+Ttj8D+H3gjcDtEXFN2jB3cp+CmVlx3epT+GPgDOB6SQcBm4FhsstOXwUuj4jVnQZbRqnq7YvstyzrArj3oGTqXvNf9+Mrog/nYtqkEBGPA/8I/GO+7Oa+wLbJ5al1lGqu/yL7Lcu6AKnOhc1S3ddIqPvxFdGnc9Fu9RFks6MK2EvSgZIOTBRT36Wqty+y37KsC+Deg5Kpe81/3Y+viD6di7amuZD0V8CFwM/ZOXV2AIcliquvUtXbF9lvWdYFcO9BydS95r/ux1dEn85Fu58UzgEOjogXRMSh+VctEwKkq7cvst+p5v/v9boA7j0ombrX/Nf9+Iro07loNymsB7akDKRMUtXbF9lvWdYFcO9BydS95r/ux1dEn87FTH0K5+YPfwKskPRldl1P4bKEsfVNqnr7Ivsty7oA7j0ombrX/Nf9+Iro07mYqU/hwmm+NyKi53d/3KdgZlZcV/oUIuI9+c5eGxGfb/oBr+0sxMFUhv6H1199B9/58SON50ueO4/r3nxMxzGY1cqXzoWV10DsAA3BkUvhxC5cHCl5H0a79xTOb3ObTSPVegNF9tucEAC+8+NHeP3Vd3QUg1mtfOlcGPt4lhAg+3fs49n2TlRgbYlpk4KkEyR9GFgg6YpJX9cAT/YkwhopQ/9Dc0KYabvZQFp5TbHt7apAH8ZMfQoPAGPAScDKSdsfBf5HqqDqqgz9D2bWhthRbHu7KtCHMdM9he8D35f0mYjY3qOYauuAkWE2tPhF3Y3+hxT7NRtYGmqdADT09G1F7D2aXzpqsb0k2r2ncLekNU1f35J0uaR9kkZYI2Xof1jy3Hkt9zHVdrOBdOTSYtvbVYE+jHaTwleALwOvz7++CHwL+BlwTZLIaujkIxZwySmHsmBkGAELRoa55JRDu9L/0O5+r3vzMU9LAK4+Mmty4mWw+Mydnww0lD3vtProsFPhT66AvRcCyv79kytKVX00bZ9CY5D0nYhY0mqbpLUR0bOOKvcpmJkV1631FCbsIenoiLgr3/lRwB75ay2rkCTNBb4J7J7/nBsi4sKmMbsDnwKOBDYBp0XE/W3GVEjR/oCqrSFQZO2Fup+LpHXgRWrXU8WR8vhKXkPfkaLHVudzMY12k8KbgP8taQ+y6bO3Am+S9Czgkim+59fASyPisXwthm9L+kpE3DlpzJnALyLieZJeB7wPOG1WRzKNomsCVG0NgSJrL9T9XCSdg36idn3CRO06PD0xpIoj5fHVeS2DosdW53Mxg7buKUTE9/JLRIcDh0fEYRHxfyLilxHRsusiMo/lT+fkX83Xql4FfDJ/fANwnNT9aUCL9gdUbQ2BImsv1P1cJK0DL1K7niqOlMdXgRr6WSt6bHU+FzNodz2F3YFXA4uA3SZ+b88095GkIbL+hucBH5m4/DTJArIZWImIJyVtAfYBHm7az1nAWQAHHlh8bZ+idfxVq/svsvZC3c9F0jrwIrXrqeJIeXwVqKGftaLHVudzMYN2q4/+leyv+ieBX076mlZE7IiIw4FR4ChJhzQNafWp4Gm/ySLiqohYHBGL58+f32bIOxVdE6BqawgUWXuh7uci6Rz0U9Wot9qeKo6Ux1fntQyKHludz8UM2k0KoxFxWkS8PyL+YeKr3R+Sr+m8AvjjppfGgYUAknYD9ga6Pt9C0f6Aqq0hUGTthbqfi6R14EVq11PFkfL4KlBDP2tFj63O52IG7SaF70oqVHYqab6kkfzxMPAy4N6mYTcBb8wfvwa4LdqpkS2oaH9Aqn6CVC4++VDe8OIDG58MhiTe8OIDW1Yf1f1cJK0DL1K7niqOlMdXgRr6WSt6bHU+FzNot0/hh2T3Bf6drKpIZPeSp1ySU9JhZDeRh8iSz7KIuEjSRcBYRNyUl61+GjiC7BPC6yLiJ9PF4j4FM7Piut2ncELRACJiDdkv++btF0x6/DjgdRnMzEqi3ZLUn5Jd+39p/vhX7X5vVS1ftYEl772Ng97xZZa897aO1zywmlizDC4/BN49kv073Tz4RcamUjSGMhxf1fZbM+2WpF4ILAYOBj5B1nNwLbBkuu+rqso1bFlvFGloKkPzU8qGrao155Xh/aiIdv/a/69kayr8EiAiHgD2TBVUv1WuYct6o0hDUxman1I2bFWtOa8M70dFtJsUnsirggIgn96itirXsGW9UaShqQzNTykbtqrWnFeG96Mi2k0KyyT9MzAi6c3A14Gr04XVX5Vr2LLeKNLQVIbmp5QNW1VrzivD+1ER7d5o/gDZ3ERfILuvcEFEfDhlYP1UuYYt640iDU1laH5K2bBVtea8MrwfFdFuSSoR8TXgawljKY2Jm8mVmi7a0pu4IdnOdMpFxpYh3qLjUx1f1fZbQ9M2r0l6lBZzEbGzeW2vVIFNxc1rZmbFdaV5LSJqW2FkllyRBXnKomoxl2UhnLLE0QVtXz4yswKKLMhTFlWLuSy9B2WJo0tq3ZVs1jdFFuQpi6rFXJbeg7LE0SVOCmYpFFmQpyyqFnNZeg/KEkeXOCmYpVBkQZ6yqFrMZek9KEscXeKkYJZCkQV5yqJqMZel96AscXSJk4JZCkUW5CmLqsVcloVwyhJHl7S1yE6ZuE/BzKy4bi+yY9Z9VaztThVzqv6AKp5j6ysnBeuPKtZ2p4o5VX9AFc+x9Z3vKVh/VLG2O1XMqfoDqniOre+cFKw/qljbnSrmVP0BVTzH1ndOCtYfVaztThVzqv6AKp5j6zsnBeuPKtZ2p4o5VX9AFc+x9Z2TgvVHFWu7U8Wcqj+giufY+s59CmZmA6DdPgV/UjBbswwuPwTePZL9u2ZZ7/ebKgazgtynYIMtVS1/kf26n8BKxJ8UbLClquUvsl/3E1iJOCnYYEtVy19kv+4nsBJxUrDBlqqWv8h+3U9gJeKkYIMtVS1/kf26n8BKxEnBBluqWv4i+3U/gZWI+xTMzAZA3/sUJC2UdLukdZLukXROizF/KGmLpNX5lz8vV10V6+3dT5Cez1tlpOxTeBJ4W0TcLWlPYKWkr0XED5vGfSsiTkwYh/VKFevt3U+Qns9bpST7pBARD0bE3fnjR4F1wIJUP89KoIr19u4nSM/nrVJ6cqNZ0iLgCOCuFi8fI+n7kr4i6QVTfP9ZksYkjW3cuDFhpNaRKtbbu58gPZ+3SkmeFCTtAXwB+OuI2Nr08t3Ab0XE7wEfBpa32kdEXBURiyNi8fz589MGbLNXxXp79xOk5/NWKUmTgqQ5ZAnhuoi4sfn1iNgaEY/lj28G5kjaN2VMllAV6+3dT5Cez1ulpKw+EvBxYF1EtJwYXtL++TgkHZXHsylVTJZYFevt3U+Qns9bpSTrU5D0+8C3gLXAU/nmdwIHAkTERyW9FfhLskqlbcC5EfHd6fbrPgUzs+La7VNIVpIaEd8GNMOYK4ErU8VgU1izLKv82DKeXdc97oLB/qvtS+fCymsgdmSrnh25tPNVz8wqyuspDBrXjO/qS+fC2Md3Po8dO587MdgA8txHg8Y147taeU2x7WY156QwaFwzvqvYUWy7Wc05KQwa14zvSkPFtpvVnJPCoHHN+K6OXFpsu1nNOSkMGteM7+rEy2DxmTs/GWgoe+6bzDagvJ6CmdkA6HufwiBZvmoDl956Hw9s3sYBI8Ocd/zBnHxEjSaErXtfQ92Prwx8jivDSaFDy1dt4Pwb17Jte1atsmHzNs6/cS1APRJD3fsa6n58ZeBzXCm+p9ChS2+9r5EQJmzbvoNLb72vTxF1Wd37Gup+fGXgc1wpTgodemDztkLbK6fufQ11P74y8DmuFCeFDh0wMlxoe+XUva+h7sdXBj7HleKk0KHzjj+Y4Tm7NjoNzxnivOMP7lNEXVb3voa6H18Z+BxXim80d2jiZnJtq48mbgTWtXKk7sdXBj7HleI+BTOzAdBun4IvH5nV2ZplcPkh8O6R7N81y6qxb+sbXz4yq6uU/QHuPagtf1Iwq6uU/QHuPagtJwWzukrZH+Deg9pyUjCrq5T9Ae49qC0nBbO6Stkf4N6D2nJSMKurlGtneF2O2nKfgpnZAHCfgpmZFeakYGZmDU4KZmbW4KRgZmYNTgpmZtbgpGBmZg1OCmZm1uCkYGZmDcmSgqSFkm6XtE7SPZLOaTFGkq6Q9CNJayS9MFU81gHPm282MFKup/Ak8LaIuFvSnsBKSV+LiB9OGnMC8Pz862jgn/J/rSw8b77ZQEn2SSEiHoyIu/PHjwLrgOaFi18FfCoydwIjkp6TKiabBc+bbzZQenJPQdIi4AjgrqaXFgDrJz0f5+mJA0lnSRqTNLZx48ZUYVornjffbKAkTwqS9gC+APx1RGxtfrnFtzxthr6IuCoiFkfE4vnz56cI06biefPNBkrSpCBpDllCuC4ibmwxZBxYOOn5KPBAypisIM+bbzZQUlYfCfg4sC4iLpti2E3An+VVSC8GtkTEg6lislnwvPlmAyVl9dES4E+BtZJW59veCRwIEBEfBW4GXgH8CPgV8OcJ47HZOuxUJwGzAZEsKUTEt2l9z2DymADekioGMzMrxh3NZmbW4KRgZmYNTgpmZtbgpGBmZg1OCmZm1uCkYGZmDU4KZmbWoKxVoDokbQR+2u84prAv8HC/g0jIx1dddT428PG147ciYsbJ4yqXFMpM0lhELO53HKn4+KqrzscGPr5u8uUjMzNrcFIwM7MGJ4XuuqrfASTm46uuOh8b+Pi6xvcUzMyswZ8UzMyswUnBzMwanBRmQdKQpFWSvtTitaWSNkpanX+9qR8xdkLS/ZLW5vGPtXhdkq6Q9CNJayS9sB9xzkYbx/aHkrZMev8qte6opBFJN0i6V9I6Scc0vV7Z9w7aOr7Kvn+SDp4U92pJWyX9ddOY5O9fypXX6uwcYB2w1xSvfy4i3trDeFJ4SURM1SxzAvD8/Oto4J/yf6tiumMD+FZEnNizaLrrQ8AtEfEaSc8Antn0etXfu5mODyr6/kXEfcDhkP3hCWwA/qVpWPL3z58UCpI0CrwS+Fi/Y+mjVwGfisydwIik5/Q7qEEnaS/gWLK10YmIJyJic9Owyr53bR5fXRwH/DgimmdvSP7+OSkU90Hg7cBT04x5df7R7gZJC3sUVzcF8FVJKyWd1eL1BcD6Sc/H821VMNOxARwj6fuSviLpBb0MrkO/DWwEPpFf3vyYpGc1janye9fO8UF137/JXgdc32J78vfPSaEASScCD0XEymmGfRFYFBGHAV8HPtmT4LprSUS8kOyj6lskHdv0equ1t6tS2zzTsd1NNkfM7wEfBpb3OsAO7Aa8EPiniDgC+CXwjqYxVX7v2jm+Kr9/AOSXxU4CPt/q5Rbbuvr+OSkUswQ4SdL9wGeBl0q6dvKAiNgUEb/On14NHNnbEDsXEQ/k/z5Edk3zqKYh48DkT0CjwAO9ia4zMx1bRGyNiMfyxzcDcyTt2/NAZ2ccGI+Iu/LnN5D9Em0eU8n3jjaOr+Lv34QTgLsj4uctXkv+/jkpFBAR50fEaEQsIvt4d1tEvGHymKbreyeR3ZCuDEnPkrTnxGPg5cAPmobdBPxZXgnxYmBLRDzY41ALa+fYJO0vSfnjo8j+H9nU61hnIyJ+BqyXdHC+6Tjgh03DKvneQXvHV+X3b5LTaX3pCHrw/rn6qAskXQSMRcRNwNmSTgKeBB4BlvYztln4TeBf8v+vdgM+ExG3SPoLgIj4KHAz8ArgR8CvgD/vU6xFtXNsrwH+UtKTwDbgdVGttv+/Aq7LL0H8BPjzmrx3E2Y6vkq/f5KeCfwR8N8nbevp++dpLszMrMGXj8zMrME3GTqVAAACwElEQVRJwczMGpwUzMyswUnBzMwanBTMzKzBScFqT9JjCfZ5uKRXTHr+bkl/M8XYYUn/lk9y1snPfIakb0pyKbkl46RgNjuHk9WLt+MM4MaI2NHJD4yIJ4BvAKd1sh+z6Tgp2ECRdJ6k7+UTFr4n37Yon5v/akn3SPqqpOH8tRflY++QdKmkH+SNUxcBp+Xz3k/8kv5dSSsk/UTS2ZN+7OuBf50Uw9uVrenwfUnvzbetkHR5/klgXf5zb5T0/yRdPGlfy/P9mSXhpGADQ9LLyeahP4rsL/0jJ02I93zgIxHxAmAz8Op8+yeAv4iIY4Ad0PiL/QKydTMOj4jP5WN/Bzg+3/+FkubkCeS3I+L+PIYTgJOBo/NJ294/KcQnIuJY4KNkSeQtwCHAUkn75GN+ALyoW+fErJmTgg2Sl+dfq8hm0/wdsmQA8O8RsTp/vBJYJGkE2DMivptv/8wM+/9yRPw6X8DnIbJpNfYlSzITXgZ8IiJ+BRARj0x67ab837XAPRHxYD654k/IJ0HLL0E9MTGHk1m3+YaVDRIBl0TEP++yUVoE/HrSph3AMK2nKZ5O8z52A7YAc5timGpumYnvf6ppX0+x6/+ruwOPF4zNrC3+pGCD5FbgDEl7AEhaIGm/qQZHxC+AR/PZKCGbGXfCo8CMf63n+xiSNJEYvprH8Mw8hnlFDiC/jLQxIrYX+T6zdjkp2MCIiK+SXQK6Q9Jasvn4Z/rFfiZwlaQ7yP7K35Jvv53sxvLkG81T+Srw+3kMt5BdJhqTtBpoWcY6jZeQzZRploRnSTWbhqQ9JhZtkfQO4DkRcU7BfRwBnBsRf9qFeG4Ezs8XeTfrOt9TMJveKyWdT/b/yk+ZxfoYEbFK0u2ShjrpVcgrmZY7IVhK/qRgZmYNvqdgZmYNTgpmZtbgpGBmZg1OCmZm1uCkYGZmDf8fULVVPLpKFrwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 展示数据\n",
    "visualization = True\n",
    "X,y = gen_data(data_size=data_size,visualization=visualization)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "参数说明：\n",
    "    ```\n",
    "    self.X:             训练数据\n",
    "    self.y:             标签\n",
    "    self.m,self.n:      m,n分别代表训练数据的行数和列数\n",
    "    self.max_iter；     最大迭代次数\n",
    "    self.C:             KKT条件\n",
    "    self.epsilon:       松弛变量\n",
    "    self.fun_kernel:   核函数选择：可选参数: linear、poly、gaussian\n",
    "    self.alpha:         核心参数\n",
    "    self.E:             核心参数\n",
    "    self.sigma:         使用高斯核时，高斯函数的参数\n",
    "    self.b:             初始化 = 0.0\n",
    "    self.w:             初始化 = 0.0\n",
    "    ```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SVM:\n",
    "    def __init__(self, max_iter=100, C=1.0, epsilon=0.001, fun_kernel='linear'):\n",
    "        self.max_iter = max_iter\n",
    "        self.C = C\n",
    "        self.fun_kernel=fun_kernel\n",
    "        # 松弛变量\n",
    "        self.epsilon = epsilon\n",
    "        # 使用高斯核时的参数sigma\n",
    "        self.sigma = 10.0\n",
    "    \n",
    "    def _init_param(self, X, y):\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "        self.b = 0.0\n",
    "        self.w = 0.0\n",
    "        # m,n分别是训练数据的行数和列数\n",
    "        self.m,self.n = X.shape\n",
    "\n",
    "        self.alpha = np.zeros(self.m)\n",
    "        self.E = [self.E_i(i) for i in range(self.m)]\n",
    "        \n",
    "        self.support_vectors = []\n",
    "    \n",
    "    # 核函数计算\n",
    "    def kernel(self, x1, x2):\n",
    "        if self.fun_kernel == 'linear':\n",
    "            return np.dot(x1, x2)\n",
    "        elif self.fun_kernel == 'poly':\n",
    "            return (np.dot(x1, x2)+1)**2\n",
    "        elif self.fun_kernel == 'gaussian':\n",
    "            return np.exp(-1 * np.matmul(x1-x2,(x1-x2).T) / (2* (self.sigma **2)))\n",
    "    \n",
    "    # 外层循环时求 g(x)\n",
    "    def gx(self, i):\n",
    "        g = self.b\n",
    "        for j in range(self.m):\n",
    "            g += self.alpha[j] * self.y[j] * self.kernel(self.X[i], self.X[j])\n",
    "        return g\n",
    "    \n",
    "    # 外层循环时求 E_i\n",
    "    def E_i(self, i):\n",
    "        return self.gx(i) - self.y[i]\n",
    "    \n",
    "    # 判断第i个alpha是否满足KKT条件，使用加了松弛变量的条件\n",
    "    def KKT(self, i):\n",
    "        y_gx = self.y[i] * self.gx(i) \n",
    "        if (abs(self.alpha[i])<self.epsilon) and (y_gx >= 1):\n",
    "            return True\n",
    "        elif (abs(self.alpha[i]-self.C) <self.epsilon) and (y_gx <= 1):\n",
    "            return True\n",
    "        elif (self.alpha[i]>-self.epsilon) and (self.alpha[i] < (self.C+self.epsilon)) and (abs(y_gx-1) <self.epsilon):\n",
    "            return True\n",
    "        else:\n",
    "            return False\n",
    "\n",
    "    # 截断公式\n",
    "    def cut_off(self, alpha, i1, i2):\n",
    "        # 计算边界\n",
    "        L=H=0\n",
    "        if self.y[i1] == self.y[i2]:\n",
    "            L = max(0, self.alpha[i1]+self.alpha[i2]-self.C)\n",
    "            H = min(self.C, self.alpha[i2]+self.alpha[i1])\n",
    "        else:\n",
    "            L = max(0, self.alpha[i2]-self.alpha[i2])\n",
    "            H = min(self.C, self.C+self.alpha[i2]-self.alpha[i1])\n",
    "        # 截断\n",
    "        if alpha > H:\n",
    "            return H \n",
    "        elif alpha <L:\n",
    "            return L\n",
    "        else:\n",
    "            return alpha\n",
    "    \n",
    "    # 两重循环选择变量alpha\n",
    "    def select_alpha(self):\n",
    "        # 外层循环首先看0<a<C的点是否满足KKT条件\n",
    "        idx_list = [i for i in range(self.m) if 0<self.alpha[i]<self.C]\n",
    "        # 否则遍历整个数据集\n",
    "        non_satisfy_list = [i for i in range(self.m) if i not in idx_list]\n",
    "        idx_list.extend(non_satisfy_list)\n",
    "\n",
    "        for i in idx_list:\n",
    "            if self.KKT(i):\n",
    "                continue\n",
    "            E1 = self.E[i]\n",
    "            if E1 >= 0:\n",
    "                j = np.argmin(self.E)\n",
    "            else:\n",
    "                j = np.argmax(self.E)\n",
    "            return i,j\n",
    "\n",
    "    # 计算阈值参数b\n",
    "    def cal_b(self, alpha1_new, alpha2_new, K11,K22,K12,K21, E1, E2, i1,i2):\n",
    "        b1_new = -E1 - self.y[i1]*K11*(alpha1_new-self.alpha[i1]) - self.y[i2]*K21*(alpha2_new-self.alpha[i2]) + self.b\n",
    "        b2_new = -E2 - self.y[i1]*K12*(alpha1_new-self.alpha[i1]) - self.y[i2]*K22*(alpha2_new-self.alpha[i2]) + self.b\n",
    "        if 0 < alpha1_new < self.C:\n",
    "            return b1_new\n",
    "        elif 0 < alpha2_new < self.C:\n",
    "            return b2_new\n",
    "        else:\n",
    "            return (b2_new+b1_new)/2\n",
    "    \n",
    "    # 训练过程\n",
    "    def fit(self, X, y):\n",
    "        self._init_param(X, y)\n",
    "        # train\n",
    "        for k in range(self.max_iter):\n",
    "            # 选择两个alpha\n",
    "            i1,i2 = self.select_alpha()\n",
    "            E1 = self.E[i1]\n",
    "            E2 = self.E[i2]\n",
    "            K11 = self.kernel(self.X[i1],self.X[i1])\n",
    "            K12 = self.kernel(self.X[i1],self.X[i2]) \n",
    "            K21 = self.kernel(self.X[i2],self.X[i1]) \n",
    "            K22 = self.kernel(self.X[i2],self.X[i2])\n",
    "            eta = K11 + K22 - 2*K12\n",
    "            if eta <= 0: continue\n",
    "            # 更新alpha2\n",
    "            alpha2_new_unc = self.alpha[i2] + self.y[i2] * (E1-E2) / eta\n",
    "            alpha2_new = self.cut_off(alpha2_new_unc, i1, i2) \n",
    "\n",
    "            #更新alpha1\n",
    "            alpha1_new = self.alpha[i1] + self.y[i1] * self.y[i2] * (self.alpha[i2] - alpha2_new)\n",
    "\n",
    "            # 更新b和E\n",
    "            b_new = self.cal_b(alpha1_new, alpha2_new, K11,K22,K12,K21, E1, E2, i1,i2)\n",
    "            \n",
    "            # 更新所有参数\n",
    "            self.alpha[i1] = alpha1_new\n",
    "            self.alpha[i2] = alpha2_new\n",
    "            self.b = b_new\n",
    "            # 更新差值E\n",
    "            self.E[i1] = self.E_i(i1)\n",
    "            self.E[i2] = self.E_i(i2)\n",
    "        for i in range(self.m):\n",
    "            if self.alpha[i] > 0:\n",
    "                self.support_vectors.append(i)\n",
    "        return \"train done!\"\n",
    "\n",
    "    def predict(self, data):\n",
    "        gx = self.b\n",
    "        for i in range(self.m):\n",
    "            gx += self.alpha[i] * self.y[i] * self.kernel(data, self.X[i])\n",
    "        if gx>0:\n",
    "            fx = 1\n",
    "        else:\n",
    "            fx = -1\n",
    "        return fx\n",
    "    # 验证SVM模型在测试数据集上的准确率\n",
    "    def score(self, X_test, y_test):\n",
    "        count=0\n",
    "        for i in range(len(X_test)):\n",
    "            pred = self.predict(X_test[i])\n",
    "            if pred == y_test[i]:\n",
    "                count += 1 \n",
    "        accuracy = count/len(X_test)\n",
    "        return accuracy\n",
    "    # 返回权重w\n",
    "    def weights(self):\n",
    "        y_x = self.y.reshape(-1,1)*self.X\n",
    "        self.w = np.dot(y_x.T, self.alpha)\n",
    "        return self.w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "#一些超参数\n",
    "C = 1.0\n",
    "epsilon = 0.001\n",
    "max_iter = 70\n",
    "data_size = 100\n",
    "visualization = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "测试线性核在iris数据集结果，并可视化分类界面等信息"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "支持向量数为 6\n",
      "svm预测准确率： 1.0\n",
      "w权重分别表示为 1.6999999999999984 -0.9000000000000004 -6.414981825659044\n",
      "可视化分类界面如下\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEWCAYAAABliCz2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4VGX2wPHvSSMh9F5CKIL0ABopFpS1K5a1t12xrrv6c90F67qo2FZRURfXVeyra6GIXVBsq2IB6b1DqKEklPTk/P6YCYYwITfJ3LlTzud58pDcubn3vAycvDlz5n1FVTHGGBP94rwOwBhjTGhYwjfGmBhhCd8YY2KEJXxjjIkRlvCNMSZGWMI3xpgYYQnfRAURWSQiJxzi8btE5IUg3OcVEXmgrtcxxgsJXgdgTDCoau9qHn8oVLGUE5GvgNdVtc4/aIwJBpvhm6gnIjaxMQZL+CZKiMhaETnJ//m9IjJJRF4Xkd3ACP+x1/2PJ/sf2yEiOSLys4i0ruK6A0TkFxHZIyJvA8kVHmsqIh+KSLaI7PJ/nuZ/7EHgOGC8iOwVkfH+40+JyAYR2S0is0XkOHf/Zoz5lSV8E63OASYBTYA3Kj12JdAY6AA0B24A8itfQESSgKnAf4BmwETg/AqnxAEvAx2BdP81xgOo6t+A/wE3qWoDVb3J/z0/A/391/svMFFEkjEmBCzhm2g1U1WnqmqZqlZO5sX4En1XVS1V1dmqujvANQYDicCTqlqsqpPwJWwAVHWHqk5W1TxV3QM8CBx/qKBU9XX/95Wo6uNAPaB7HcZpjGOW8E202nCIx/4DTAPeEpFNIvKoiCQGOK8dsFEPXGFwXfknIlJfRJ4TkXX+0tE3QBMRia/qxiIyUkSWiEiuiOTg+02jRU0GZkxtWcI30arKZWD9s/X7VLUXcDQwHPh9gFM3A+1FRCocS6/w+Uh8s/NBqtoIGOo/Xn7+ATH46/W3AxcBTVW1CZBb4XxjXGUJ38QcERkmIn39M/Hd+Eo8pQFOnQmUADeLSIKInAcMrPB4Q3x1+xwRaQbcU+n7twJdKp1fAmQDCSIyGmgUjDEZ44QlfBOL2uB7QXc3sAT4Gni98kmqWgScB4wAdgEXA1MqnPIkkAJsB34APq10iaeAC/wdPE/jKyN9AizHVxoq4NClJ2OCSmwDFGOMiQ02wzfGmBhhCd8YY2KEJXxjjIkRlvCNMSZGhNWiUi1atNBOnTp5HYYxxkSM2bNnb1fVlk7ODauE36lTJ2bNmuV1GMYYEzFEZF31Z/lYSccYY2KEJXxjjIkRriZ8EWniX5d8qX/BqCFu3s8YY0zV3K7hPwV8qqoX+NcWr+/y/YwxUay4uJisrCwKCgq8DiXkkpOTSUtLIzEx0MKuzriW8EWkfPXAEbB/XZIit+5njIl+WVlZNGzYkE6dOnHgIqbRTVXZsWMHWVlZdO7cudbXcbOk0wXfqoAvi8gcEXlBRFIrnyQi14vILBGZlZ2d7WI4xphIV1BQQPPmzWMq2QOICM2bN6/zbzZuJvwE4AjgWVUdAOwD7qh8kqo+r6qZqprZsqWjVlJjTAyLtWRfLhjjdjPhZwFZqvqj/+tJ+H4AGGNM1CopLWNTTj4lpWVeh3IQ1xK+qm4BNohI+X6dJwKL3bqfMcZ4bU9BMSu27WXH3iL2FQXaUyewtWvX0qdPHwDmzp3Lxx9/7Ep8bvfh/x/whojMB/oDD7l8P2OMCRlVpaysjNKyMrJ25bFm+z7iRTisVSqNU2rXTROxCV9V5/rr8xmqeq6q7nLzfsYY47a1a9fSs2dP/vSnP3HEEUcw4aWXOXLgYE4eOoS/3XQ1bepD/aQE7rjjDnr16kVGRgajRo0CYMSIEUyaNGn/tRo0aHDAtYuKihg9ejRvv/02/fv35+233w5q7GG1lo4xxjh13weLWLxpd1Cv2atdI+45q3e15y1btowXXnyJP/71TkZcfjEvv/M+h7dvwfgnH+fJJ8dx00038e6777J06VJEhJycHEf3T0pKYsyYMcyaNYvx48fXdTgHsYRvjDE1lJ7ekeZdevPZpx+zZuVyrjrvNMA3Qx8yZAiNGjUiOTmZa6+9ljPPPJPhw4d7HLGPJXxjTERyMhMPtrIyZevuAhLqJQPQulEyp55yMm+++eZB5/7000/MmDGDt956i/Hjx/PFF1+QkJBAWZmve0dVKSoK7XtRbfE0Y4xxYF9hCSu27WXnviIS4uLo1qohw447hu+++46VK1cCkJeXx/Lly9m7dy+5ubmcccYZPPnkk8ydOxfwLQE/e/ZsAN577z2Ki4sPuk/Dhg3Zs2ePK2OwhG+MMYdQpsrm3HxWZ+9FVUlvlkJCvBAfJ7Rs2ZJXXnmFSy+9lIyMDAYPHszSpUvZs2cPw4cPJyMjg+OPP55x48YBcN111/H1118zcOBAfvzxR1JTD1p8gGHDhrF48WJXXrQVVQ3qBesiMzNTbQMUY0xVlixZQs+ePUN2v7yiErJ25VNQXEqz1CTaNk4mPs67eXKg8YvIbFXNdPL9VsM3xphKylTZtqeQ7N2FJMQLnVqk0ii59qtUhgtL+MYYU0F+USlZu/LILy6laX3frD4hPjqq35bwjTEGX9dM9p5Ctu4pJF6Ejs1r/27ZcGUJ3xgT8wqKfbP6vKJSGqck0r5JStTM6iuyhG+MiVmqyva9hWzZXUi8QHqz+jSpn+R1WK6xhG+MiUmFxaVk7cpnX1EJjZITad80hcQonNVXFN2jM8aYSspn9Su27aWgpJQOzerTsXn9Gif7pUuXMmTIEOrVq8djjz3mUrTBZTN8Y0zMKCopZcOufPYVltAwOZG0JikkJtRu3tusWTOefvpppk6dGuQo3WMJ3xgTtabO2cjYacvYlJNP68bJXD4ond90b0Va0xSa1k+q07aBrVq1olWrVnz00UdBjNhdlvCNMVFp6pyN3DllAfnFvp2ntuQWMP6LlbRtlEzv9o09js4bVsM3xkSlsdOW7k/25QpLyhj3+QqPIvKeJXxjTNQpLi1jU05BwMc25eTX+rrPPPMM/fv3p3///mzatKnW1/GKlXSMMVFDVcnNL2ZjTj4tGtYje0/hQee0a5JS6+vfeOON3HjjjXUJ0VM2wzfGRIXi0jLW78xj/c486iXEc9up3UlJjD/gnJTEeG49tXtQ7rdlyxbS0tJ44okneOCBB0hLS2P37uBuuRhsNsM3xkS83LwiNuYUUKpKm8bJtGxQj66tGpAYH7e/S6ddkxRuPbU75w5oH5R7tmnThqysrKBcK1Qs4RtjIlaJv1afk19ESmI8XZqlklxhVn/ugPZBS/DRwBK+MSYi7c4vJisnn9JSpXWjZFo2rEdcHfrqY4ElfGNMRCkpK2NzTgG78opIToync/P6pCRZKnPC/paMMRGjoLiUFVv3UlKqtGqYTKtGNquvCevSMcaEvb2FJdw5ZQHb9xYRJ8JhrVJp0zjZkn0NWcI3xoS171dt59Rx3/DWz+tpmJxAt1YNqG8lnFpxNeGLyFoRWSAic0Vklpv3MsZEl7yiEu55byGXTfiRpIQ4Jt0whMYpicTFeTurv/rqq2nVqhV9+vTxNI7aCMUMf5iq9lfVzBDcyxgTBWat3ckZT/2PV2euY8TRnfj45uM4smMzr8MCYMSIEXz66adeh1Er9nuRMSZsFBSX8vj0Zbzw7RrSmqbw1vWDGdylee0vOP8dmDEGcrOgcRqcOBoyLqpTjEOHDmXt2rV1uoZX3E74CkwXEQWeU9XnK58gItcD1wOkp6e7HI4xJlzN3ZDDyHfmsip7H5cPSueuM3qSWq8OKWr+O/DBzVDsXywtd4Pva6hz0o9Ubif8Y1R1k4i0Aj4TkaWq+k3FE/w/BJ4HyMzMVJfjMcaEmcKSUp6esYJnv1pF60bJ/OeagRzXrWXdLzxjzK/Jvlxxvu+4JfzgU9VN/j+3ici7wEDgm0N/lzEmVizcmMuoifNYumUPFxyZxuizetEoOTE4F8+tYp2bqo7HANcSvoikAnGqusf/+SnAGLfuZ4yJHMWlZTzz5UrGf7GSpqlJvHhlJif2bB3cmzRO85VxAh2PUW526bQGvhWRecBPwEeqGpkvbRtjgmbZlj389l/f8eTnKzgzoy2f/WVo8JM9+F6gTay09n1iiu94HVx66aUMGTKEZcuWkZaWxosvvlin64WSazN8VV0N9HPr+saYyFJSWsbz/1vNk5+toGFyAv++4ghO69PWvRuW1+mD3KXz5ptvBiE4b1hbpjHGdauy9zJq4jzmrM/hjL5tuP+cPjRvUM/9G2dcFLMv0AZiCd8Y45qyMuWl79YwdtoyUpLiefrSAZyV0RaxNXA8YQnfGOOKdTv2cevE+fy0dicn9WzFQ+f1pVXD5DpfV1Vj8geGat271i3hG2OCqqxMef3HdTz88VIS4oXHLuzH+Ue0D0qSTk5OZseOHTRv3jymkr6qsmPHDpKT6/YD0xK+MSZosnblcduk+Xy/agdDD2/JI+f3pW3jlOq/0aG0tDSysrLIzs4O2jUjRXJyMmlpdWsptYRvjKkzVeXtnzfwwEdLUFUePq8vlxzVIeiz8MTERDp37hzUa8YSS/jGmDrZnJvPHZMX8PXybIZ0ac6jF2TQoVl9r8MyAVjCN8bUiqoy5ZeN3PvBIkpKlTHn9OaKQR09X6/eVM0SvjGmxrbtKeCuKQv5fMlWjurUlLEX9KNTi1SvwzLVsIRvjHFMVflg/mZGv7eQvKJS7j6zJ1cd05l4m9VHBEv4xhhHduwt5O/vLeTjBVvo16EJj1/Yj66tGngdlqkBS/jGmGp9unALf3t3AXsKSrjttO5cf1wXEuJDsUOqCSZL+MaYKuXkFXHv+4uYOncTvds14r/X9ad7m4Zeh2VqyRK+MSEydc5Gxk5bxqacfNo1SeHWU7tz7oD2XodVpS+WbuWOyQvYua+IW07qxo3DupJos/qIZgnfmBCYOmcjd05ZQH5xKQAbc/K5c8oCgLBL+rsLirn/g8VMnJ1F99YNeWnEUfRp39jrsEwQWMI3JgTGTlu2P9mXyy8uZey0ZWGV8L9Zns3tk+ezdXcBfzrhMP58UjfqJcR7HZYJEkv4xoTAppz8Gh0Ptb2FJTz08RL+++N6urRMZfIfj2ZAelOvwzJBZgnfmBBo1ySFjQGSe7smwVtYrLZmrtrBrZPmsTEnn2uP7cyoU7uTnGiz+mhkr8AYEwK3ntqdlEpJNCUxnltP7e5RRJBfVMq97y/i0gk/EB8nvPOHIdw9vJcl+yhmM3xjQqC8Th8uXTqz1+1k1MT5rNm+jyuHdOT203tQP8nSQbSzZ9iYKgS7jfLcAe09f4G2oLiUcZ8tZ8L/VtO2cQr/vW4QRx/WwtOYTOhYwjcmgEhqo3RqflYOI9+Zx4pte7l0YDp/O7MnDepZCogl9mwbE0CktFE6UVRSxj+/WMG/vlpFywb1ePXqgRx/eEuvwzIesIRvwp4X71AN9zZKpxZv2s3IifNYsnk35x+RxuizetE4JdHrsIxHLOGbsOZVaSWc2yidKC4t49mvVvH0jBU0qZ/EhN9ncnKv1l6HZTxmbZkmrB2qtOKmcGyjdGr51j2c96/veeKz5ZzRty2f/WWoJXsD2AzfhDmvSivh1kbpRGmZMuF/q3li+nIaJCfw7OVHcHrftl6HZcKI6wlfROKBWcBGVR3u9v1MdPGytOJlG2VNX7dYnb2XkRPnMWd9Dqf1bsMDv+1Diwb1QhixiQShKOn8GVgSgvuYKBTJpZXaKn/dYmNOPsqvr1tMnbPxoHPLypQXv13D6U/9j9XZ+3jqkv48e8URluxNQK4mfBFJA84EXnDzPiZ6nTugPQ+f15f2TVIQoH2TFB4+r29Yl1bqyunrFut35HHJhB+4/8PFHNO1BdP/MpRz+rdHxPaXNYG5XdJ5ErgNqHKLHBG5HrgeID093eVwTCQKh3eohlJ1r1uoKq//uJ6HP15CvAiPXpDBhUemWaI31XIt4YvIcGCbqs4WkROqOk9VnweeB8jMzFS34jGm3N1TF/DmjxsoVSVehEsHdeCBc/t6HdZ+h3rdYmNOPrdPms+3K7dzXLcWPHJ+RsS0ihrvuVnSOQY4W0TWAm8BvxGR1128nzHVunvqAl7/YT2l6ptblKry+g/ruXvqAo8j+1Wg1y2SE+I4rlsLThv3Db+s38WDv+3Da1cPtGRvasS1hK+qd6pqmqp2Ai4BvlDVK9y6nzFOvPnjhhod90Ll1y3aNEqmS8sGvPXzBnq1a8S0W4Zy+aCOVsIxNWZ9+MYzl0+YyXerdu7/+pjDmvHGdUNcvWf5zN7p8WCqSavluQPac07/drw3dxP3vL+I1dv3cs9ZvbhySCfi4izRm9qp0QxfRFL9ffU1oqpfWQ++qahysgf4btVOLp8w09X7xlcxK67qeLDUpNUSIHtPITe8Pptb3p5L11YN+Pjm47jqmM6W7E2dHDLhi0iciFwmIh+JyDZgKbBZRBaJyFgR6RaaME20qZzsqzseLJcO6lCj48FSkyUiPpq/mVPGfc2Xy7K564wevPOHIXRp2cDV+ExsqK6k8yXwOXAnsFBVywBEpBkwDPiHiLyrqvZirIkI5d04oe7ScbJExM59RYx+byEfzt9Mv7TGPHZhP7q1rrKj2Zgaqy7hn6SqxZUPqupOYDIwWURsrVUTFpzWyDM7NuPLpdlsysmnTeNkMjs2cz226paImL5oC3e9u5Dc/CJGnXI4Nxx/GAnxtrahCa5DJvyKyV5EmgIdKn6Pqv4S6AeCMdXp1iqVFdv2BTxeG06XUfZqueVbT+1+wH3Bt0TEjcMO469vz2XKnI30bNuI164eSK92jVyLw8Q2R106InI/MAJYBZS3MyjwG3fCMtEur6isRser43SHKq92sgq0+ubwjLY8NWMF2/cWcfOJ3bhpWFeSEmxWb9zjtC3zIuAwVS1yMxhTPS92f6opJzEGe9ljp9fzcier8iUi9hQU8+BHS3jum9Uc3roBL/z+KPqmNXb9/sY4TfgLgSbANhdjMdWIhI21ncYY7GWPnV7P652svl+5nVsnzWdzbj43HH8Yfzm5G/USatzpbEytOP398WFgjohME5H3yz/cDMwczKvdn2rCaYzBXvbY6fW8Wm45r6iE0e8t5LIXfqReQhyT/ng0d5zew5K9CSmnM/xXgUeABUDtiqymziJhY22nMQZ7Rymn1/NiJ6uf1uxk1MR5bNiVxzXHdmbUKd1JSbJEb0LPacLfrqpPuxqJqZbX5QgnahLjrHU72ZJbgAJbcguYtW5nwMTr9HULp8soOz2vrq+XFPh/s3npuzV0aFqft64bzKAuzR1/vzHB5jThzxaRh4H3gcLyg6r6iytRmYCqau0Lp92fOjUPnPA7NT8w4ZevWlmufNVK4IA3QXn1ukVd7ztn/S5GTpzH6ux9/G5wR+44vQep9WzpKuMtp/8CB/j/HFzhmLVlhlgkbKz9w+pdjo4fatXKignfqzbK2t63sKSUJz9fwXNfr6Jt4xTeuHYQx3Rt4VqcxtSEo4SvqsPcDsQ4E+67PzldjdLpeV69blGb+y7IymXkxLks37qXS47qwN/O7EnDZHsjugkfTt949RDwqKrm+L9uCoxU1bvdDM4czI0+/GDuABUvEjCZV16N0ul5brxu4WS8NblvUUkZ479cyTNfrqRFgyRevuoohnVvVev4XDH/HZgxBnKzoHEanDgaMi7yOioTYk7bMk8vT/YAqroLOMOdkExVarrErhPB3gGqS8v6jo4P7tI04HmVjw/r0TLgeVUdr47T8Tpt31yyeTfnPvMdT89Ywdn92jH9luPDM9l/cDPkbgDU9+cHN/uOm5jiNOHHi0i98i9EJAWod4jzjQvc6MMP9g5Qq7PzHB1fuyNwaaTy8S+XZgc8r6rj1XE63sq7TrVvksLD5/Xd/9tUSWkZ479Ywdnjv2XbngKe/92RjLu4P43rh2EJZ8YYKK70912c7ztuYorTF21fB2aIyMv4Xqy9Gl9vvgkhN+rZwd4BKti1+WCPuSbjrer1kpXb9jDynXnMy8rlzIy23H9OH5qlJtUqnpDIzarZ8WCyUlJYcfqi7aMiMh84CRDgflWd5mpk5iBu1LOd1tKdEn5dXa/y8Yqa1E9kV97BC602qTRDDvaY6zLe0jLlxW9X89j05aQmxTP+sgEMz2hXqzhCqnGav5wT4LibyktJ5b9dlJeSwJK+R6rb8Wr//wJV/VRVR6nqyIrJXmwn5ZBxY1mAYO8AVb+Kd5BWPl7VLxCVjwd7zLUd75rt+7j4uZk89PFSjj+8JdP/cnxkJHvwzaoTK/2ATEzxHXeTlZLCTrU7XonIZOA9Vd3/LhkRSQKOBa7EtyvWK65FaPZzow8/2DtA5RWVOjqemx94G4XKx4M95pqOt6xMeW3mWv7x6VKS4uMYd3E/zu3fnoia55TPpkNdWvGylGQCqi7hn4avXv+miHQGcoAUfL8ZTAfGqepcd0OMbF4tZ1yTVssHzu0btC3+3Fi1MtjvPXA63g0787ht0nxmrt7BCd1b8o/zMmjTODlocQSF0xp5xkWhL6PUpJQUi7V+D8Zc3Y5XBcC/gH/5tzJsAeRXbNE0VQv2sgBOr+d02QI3DOvR8oB7VzxeUTgvE6GqvPnTBh78aDEiwiPn9+WizA7hN6sP9xr5iaMPjA8Cl5LCfRxu8GjMNdlepwzfa2+NRCRdRNJdiilqBLuN0un1gt1qWRNO2yira3v0yqacfH7/0k/c9e4CBqQ3ZdpfhnLxUenhl+wh/GvkGRfBWU9D4w6A+P486+mDE1q4j8MNHo3Z6Ttt/w+4B9jKr8sjK5DhUlxRwatdnWraaum07HT5hJl8t2rn/q+POawZb1w3pFYxQngtE6GqTJqdxZgPFlOqyv3n9uGKQWGa6MtFQo3cSSnJy3F8+FeY/QpoKUg8HDkChj9Rt2s6KdV4NGanM/w/A91Vtbeq9vV/WLKvRlWtg7VtKazcsljV8apaDAMdd/ru3crJHuC7VTu5fMLMA44Fe8yhsG13Ade+OotbJ82nZ7tGfPrnofxucMfwTvZQdVul2+2WwebVOD78K8x60ZfswffnrBd9x2vL6buaPRqz04S/Ach1M5BoFOyWQqetjDVpPXRaJqqc7Ks67tWOUrWhqrw3dyMnj/uGb1du5+/De/HWdYNJbx54eYiw41W7ZbB5NY7Zr9TsuBNOSzUejfmQJR0RKf9Rtxr4SkQ+4sD18Ov4u090C3ZLodNWxpq0Hga77BQJSzgDbN9byN3vLuTTRVsYkN6Exy/sR5eWDbwOq2a8arcMNq/GoYFbiKs87oTTUo1HY66uht/Q/+d6/0eS/wMCv6FyPxFJBr7Bt+ZOAjBJVe+pfajhxYt2y5q0MjptPYyEXbSC7ZMPJ/O370rYq/W4o8FnXHf08cS3PKb2F3RaB46E1kOvYgx226iTcUh84OQuddh+siatqB60ylbXlnkfgIhcqKoTKz4mIhdWc+1C4Dequtff0vmtiHyiqj/UKeIw4LQ9MthtmW60Mjrdoapbq1RWbNt30HndWqUe8LVXO1Q5sWtfEfe89gnvr0umr6zm8aRnObxkI3w02VfcrM1/vvI6cLnyOjAcmPTdaMML9jWjpT3S6Tg6HQtrvj74+zsdW/t7O21F9YjTGv6dDo/tpz57/V8m+j9qtyJXmHFa9w52W6YbrYxOd6jKKwq8d33l426s6BkMny/eyilPfsPH64S/JkxkStI9HB7nf2G6Lu1wTuvAbrThBfua0dIe6XQcO1cH/v6qjjvhtBXVI9XV8E/Ht+59exGpuIl5I6CkuouLSDwwG+gKPKOqPwY453rgeoD09Mho7fdqpUcIfitjuK9uWVe5+cWM+WAxk3/JokebhrxSOJLecWsDnFjLdjindWA32vCCfc1IaPN0wuk43BqvF+9qdqi6Gf4mYBZQgC9xl3+8D5xa3cVVtVRV+wNpwEAR6RPgnOdVNVNVM1u2rN2mFqHmtPUwEloUnbZwRuKYv16ezWlPfsPUuRu5aVhX3r/pWHo3rSJB17Ydrqp6b+XjbrTh1eSa89+BcX3g3ia+PwNtfhLs69VEMK/ndBzR0tZaA4dM+Ko6T1VfBbqq6qsVPqb4d71yxL8Uw1f41uaJeE5bD4O9W5MbnLZwOh1zOLRl7i0s4c4pC7jypZ9IrZfAlD8ezahTu5OUEBf8driq6r2Vj3c7JfB5VR13wulYnPaGO40x2DtoBft6TscRLW2tNeB0A5RfRKTy7/65+Gb/D6jqjsrfICItgWJVzfHvkHUS8Eidog0TTlsPg71bkxuctnA6HbPXbZkzV+3g1knz2JiTzx+GduEvJx9OcsUfQMFuh3NaB14xPfB5VR13wulYDlXTrniu0xidXs+pYF/P6Tiipa21Bpwm/E+AUuC//q8vwbeuTi6+pZHPCvA9bYFX/XX8OOAdVf2wTtGGESe19HCrZ1fFaQun09cPvFgyIa+ohEc/XcYr36+lU/P6TLphCEd2bBb45GDWWCOhXhzsGL187SDYyxY4/bcQCS21DjhN+MeoasVG5QUi8p2qHiMiVwT6BlWdDwyoc4QRLBZ73L0wa+1ORk2cx9odeYw4uhO3n9aDlCo2Ygk6p33XifWh+OC2VhJD8K7eYMcY7B20nF7PabtlsOOLlnZVnLdlNhCRQeVfiMhAoPxtidV268SqcKhnR7OC4lIe/GgxFz43k1JV3rp+MPee3Tt0yR6c14FLqvitrqrjwRTsGINd+3Z6Pa+WLYiWdlWcz/CvBV4SkQb4Sjm7gWtFJBV42K3gIp3X9exoNndDDiPfmcuq7H1cPiidu87oSWo9p/+cg8hpHVgDv4+hyuPBFOwYg137dno9r5YtiJZ2VZxvYv4z0FdEGgNSaQOUOvZjRbdwWgI4YlWonxY26sjTLe/l2cVJtG6UzH+uGchx3cKn66lKNX0bvxvL9lanJjEGu9fcyfW8WrbAq03gXeCopCMi9UTkMuBG4GYRGS0i0du7ZMJHhZa9hWXpnJN9A88sSuL8zsVM+8tQ75O905bCI0cE/v5Ax4O9bK8bMXrBqzbKKGrfdFrDfw84B1+9fl+FD2PcNWMMxUVFPFlyHucn05uXAAAUq0lEQVQW3c8ObciLiWMZu+9uGiUH3h8g1PE5qu8OfwIyr/l1tizxvq8DzdqDvWyvGzF6watlC8J8uYSacFr0TFPVqHjTlIksy3bByOIxLNTOnBP3HfclvkIT2Qe5YbI5SU3qu8OfcJY8g71srxsxesWrZQvCeLmEmnA6w/9eRNzd/dqYCkpKy3jmy5WcVfQgm7UZ/04cx1NJz/iSPYSmfurVcgROl2twKgaXEDCBOU34xwKzRWSZiMwXkQUiMt/NwEzsWrltLxf8eyZjpy3jpA4wPXU0p8X//OsJoaifOq17B3t5A3C+XINTUVSDNnXjtKRzuqtRGAOUlikvf7eGsdOWkZIUz9OXDuCsjLbIgsLQv8vR6dv9g728AQR/2d4YXELABOa0LXOdiBwLdFPVl/3r5ETYfnAm7FRot1yXmsGouFH8nB3PST1b8dB5fWnVMNl3nhf1Uy+XD3Cj79urGnSULEkQLRwlfBG5B8gEugMv49vM5HWgDvvCmZjmL3GUFRXweulJPLzjUhIo5LHB8Zx/TiZSxbLNIZPSFPIDbNye0vTAr914u3+09H1H0ZIE0cJpDf+3wNn4WzFVdRO/7ndrTM3NGMOGwlSuKL6T0SVXcVTcMqbXu40L1oz2PtnXhBtv94+WmnsULUkQLZzW8ItUVcuXSPYvqWBMragqb+3oxgMllwPwcMIELon/EhEg1/E2C+7KryKOysfdeLt/tNTco2hJgmjhNOG/IyLPAU1E5DrgamCCe2GZsBSEeuzm3HzumLyAr0uuZUjcIh5NeI4Ocdt/PSFcyhZOyypuvd0/Gvq+o6U0FUUclXRU9TFgEjAZXx1/tKr+083ATJip465Eqsrk2VmcMu4bflqzkzGZhbxRf9yByT6cyhbNujg7Hi3lFzfY303Ycby8oKp+BnzmYiwmnNVhV6Jtewq4a8pCPl+ylaM6NWXsBf3o1CIV5peEb9li7bfOjkdL+cUN9ncTdg6Z8EVkD1B5a0PwLZGsqtrIlahi0NQ5G8N7GeVa7kr0Qb0z+Hv+ZeSVxXP3mT256pjOxMf5X5QN57JFTZY3COdxGFPBIRO+qlonTghMnbORO6csIL/Yl0w25uRz55QFAOGT9Gu4K9GOogRGF/8fHxUMpl/cah4/vQ1dj6uiTBKOJC7w+vDitLHNWFtm+LF/vWFg7LRl+5N9ufziUsZOW+ZRRAHUYFeiaQW9ObXwUaaXZXJrwltMThxN11n3hS7WYEioYhvKqo6bg1lbZtjxYIsgU1lEbHbuoB6bm1fMPdlnMbXsWHrLGl5PfIgecf7fCiKtFa84r2bHzcGsLTPsWMIPAxGz2fkhatVfLt3G7ZPns7NsCLckTOLG+PdIlAq/tURaK160tRQGe4kDJ9eLtr/DKGAlnTAQyZud7y4o5rZJ87jqlZ9pWj+JqacVcUvKJwcm+0hsxYumlsI6ttTW+nrR9HcYJWyGHwYidbPzb1ds57ZJ89iyu4A/nXAYfz6pG/US4qFZWeS34kVTS2EdWmrrdL1o+juMEqIaqOvSG5mZmTpr1iyvwzDV2FdYwkMfL+GNH9fTpWUqj1/YjwHpTav/RuONe5tQZXf1vTneX8/UiYjMVtVMJ+faDN/UyA+rd3DrpHlk7crnuuM6M/KU7iQn1nInJhMawa6lW20+YlkN3ziSX1TKfR8s4pLnfyBOhHf+MIS/ndnLkn0kCHYt3WrzEctm+KZas9ftYtTEeazZvo8rh3Tk9tN7UD/J/ulEjGDX0q02H7Fcq+GLSAfgNaANUAY8r6pPHep7rIYfXgqKSxn3+XImfLOato1TGHtBBkd3bVH9N9ouR8aETLjU8EuAkar6i4g0xLcJ+mequtjFe5ogmZ+Vw8h35rFi214uHZjOXWf0oGFyooNvtLfTGxOuXEv4qroZ2Oz/fI+ILAHaA5bww1hRSRn//GIF//pqFS0b1OPVqwdy/OEtnV8g2C2AxpigCUkhVkQ6AQOAHwM8dj1wPUB6enoowjFVWLxpNyMnzmPJ5t2cf0Qao8/qReMUB7P6iuzt9MaELdcTvog0wLdxyi2qurvy46r6PPA8+Gr4bsdjDlZcWsazX63i6RkraFI/iQm/z+TkXq1rdzFr2TMmbLma8EUkEV+yf0NVp7h5L1M7y7fuYeQ781iwMZez+7XjvrN70zQ1qfYXPHH0gTV8sJY9Y8KEawlfRAR4EViiqk+4dR9TO6VlyoT/reaJ6ctpkJzAvy4/gjP6tq37ha1lz5iw5eYM/xjgd8ACEZnrP3aXqn7s4j2NA6uz9zJy4jzmrM/htN5teOC3fWjRoF7wbmA7QBkTltzs0vkW31aIJkyUlSkvf7+WRz9dSnJiPE9d0p+z+7XD98uYMSba2dslY8T6HXmMmjSPn9bs5Dc9WvHweX1p3SjZ67CMMSFkCT/KqSqv/7iehz9eQrwIj16QwYVHptms3pgYZAk/im3Myef2SfP5duV2juvWgkfOzwi/XbSMMSFjCT8KqSoTZ2Vx/4eLKVXlwd/24bKB6TarNybGWcKPMlt3F3DH5Pl8uSybQZ2b8diF/ejQrL7XYRljwoAl/Cihqrw3dxP3vL+IwpJS7jmrF1cO6URcnM3qjTE+lvCjQPaeQu6euoBpi7ZyRHoTHruwH11aNvA6LGNMmLGEH+E+mr+Zv7+3kL2FJdx5eg+uPa4L8TarN8YEYAk/Qu3aV8Tf31vIh/M30y+tMY9d2I9urRt6HZYxJoxZwo9Any3eyp1TFpCbX8SoUw7nhuMPIyHetic2xhyaJfwIkptXzH0fLGLKnI30bNuI164eSK92jbwOyxgTISzhR4ivlm3j9snz2b63iJt/05WbftONpASb1RtjnLOEH+b2FBTz4EdLeOvnDXRr1YAJv88kI62J12EZYyKQJfww9t3K7dw2aT6bc/O54fjDuOWkbiQnxnsdljEmQlnCD0P7Ckv4xydL+c8P6+jSIpWJNxzNkR2beh2WMSbCWcIPMz+t2cmoifPYsCuPa47tzKhTupOSZLN6Y0zdWcIPEwXFpYydtoyXvltDh6b1eeu6wQzq0tzrsIwxUcQSfhiYs34XIyfOY3X2Pn43uCN3nN6D1Hr21BhjgsuyiocKS0p58vMVPPf1Kto2TuH1awZxbLcWXodljIlSlvA9siArl5ET57J8614uOaoDfzuzJw2TE70OyxgTxSzhh1hRSRnjv1zJM1+upEWDJF6+6iiGdW/ldVjGmBhgCT+Elmzezch35rF4827OG9Cee87qTeP6Nqs3xoSGJfwQKCkt499fr+KpGStonJLI8787klN6t/E6LGNMjLGE77KV2/YwcuJ85m3I4cyMttx/Th+apSZ5HZYxJgZZwndJaZny4rereWz6clKT4hl/2QCGZ7TzOixjTAyzhO+CNdv3cevEecxat4uTe7Xmod/2pWXDel6HZYyJcZbwg6isTHlt5lr+8elSkuLjGHdxP87t3x4R23LQGOM91xK+iLwEDAe2qWoft+4TLjbszOO2SfOZuXoHJ3RvyT/Oy6BN42SvwzLGmP3cnOG/AowHXnPxHp5TVd78aQMPfrQYEeGR8/tyUWYHm9UbY8KOawlfVb8RkU5uXT8cbM7N5/bJC/hmeTbHdG3OI+dnkNa0vtdhGWNMQJ7X8EXkeuB6gPT0dI+jcUZVmTQ7izEfLqakVLn/nN5cPqgjcXE2qzfGhC/PE76qPg88D5CZmakeh1OtbbsLuHPKAmYs3cbATs0Ye2EGHZuneh2WMcZUy/OEHylUlffnbWL0e4soKC7l78N7cdXRnWxWb4yJGJbwHdi+t5C/T13IJwu3MCC9CY9d2I/DWjbwOixjjKkRN9sy3wROAFqISBZwj6q+6Nb93PLJgs3cPXUhewpKuP20Hlw/tAvxNqs3xkQgN7t0LnXr2qGwa18R97y/iPfnbaJv+8Y8flE/Dm/d0OuwjDGm1qykE8CMJVu5Y8oCdu0r4q8nH84fTziMxPg4r8Myxpg6sYRfQW5+MWM+WMzkX7Lo0aYhr1x1FL3bNfY6LGOMCQpL+H5fL8/mjsnz2bankJuGdeXmE7uRlGCzemNM9Ij5hL+3sIQHP1rMmz9toGurBky54kj6dWjidVjGGBN0MZ3wv1+1nVsnzmdTbj5/GNqFv5x8OMmJ8V6HZYwxrojJhJ9XVMIjnyzl1Znr6NS8PpNuGMKRHZt5HZYxxrgq5hL+z2t3cuvEeazdkceIoztx+2k9SEmyWb0xJvrFTMIvKC7l8enLeOHbNbRvksKb1w1myGHNvQ7LGGNCJiYS/pz1uxg1cR6rsvdx2aB07jqjJw3qxcTQjTFmv6jOeoUlpTz1+Qr+/fUqWjdK5j/XDOS4bi29DssYYzwRtQl/4cZcRk2cx9Ite7jwyDT+flYvGiUneh2WMcZ4JuoSfnFpGc98uZLxX6ykWWoSL16ZyYk9W3sdljHGeC6qEv6yLXsYOXEuCzfu5tz+7bj37N40qZ/kdVjGGBMWoiLhl5SW8dw3q3nq8xU0TE7g31ccyWl92ngdljHGhJWIT/i5ecX8/uWfmLchhzP6tuH+c/rQvEE9r8MyxpiwE/EJv1FKAp2a1+eaYztzVkZbRGxzEmOMCSTiE76I8NQlA7wOwxhjwp6t/2uMMTHCEr4xxsQIS/jGGBMjLOEbY0yMsIRvjDExwhK+McbECEv4xhgTIyzhG2NMjBBV9TqG/UQkG1hXy29vAWwPYjheiZZxQPSMJVrGAdEzlmgZB9R9LB1V1dFGH2GV8OtCRGapaqbXcdRVtIwDomcs0TIOiJ6xRMs4ILRjsZKOMcbECEv4xhgTI6Ip4T/vdQBBEi3jgOgZS7SMA6JnLNEyDgjhWKKmhm+MMebQommGb4wx5hAs4RtjTIyIuIQvIvEiMkdEPgzwWD0ReVtEVorIjyLSKfQROlfNWEaISLaIzPV/XOtFjE6IyFoRWeCPc1aAx0VEnvY/L/NF5Agv4qyOg3GcICK5FZ6T0V7E6YSINBGRSSKyVESWiMiQSo9HynNS3Tgi4jkRke4VYpwrIrtF5JZK57j+nETijld/BpYAjQI8dg2wS1W7isglwCPAxaEMroYONRaAt1X1phDGUxfDVLWqN4+cDnTzfwwCnvX/GY4ONQ6A/6nq8JBFU3tPAZ+q6gUikgTUr/R4pDwn1Y0DIuA5UdVlQH/wTfSAjcC7lU5z/TmJqBm+iKQBZwIvVHHKOcCr/s8nASdKmG5y62As0eQc4DX1+QFoIiJtvQ4qWolII2Ao8CKAqhapak6l08L+OXE4jkh0IrBKVSuvKuD6cxJRCR94ErgNKKvi8fbABgBVLQFygeahCa3GqhsLwPn+X+0miUiHEMVVGwpMF5HZInJ9gMf3Py9+Wf5j4aa6cQAMEZF5IvKJiPQOZXA10AXIBl72lwxfEJHUSudEwnPiZBwQGc9JRZcAbwY47vpzEjEJX0SGA9tUdfahTgtwLOz6Th2O5QOgk6pmAJ/z628u4egYVT0C36+kN4rI0EqPR8TzQvXj+AXfuiX9gH8CU0MdoEMJwBHAs6o6ANgH3FHpnEh4TpyMI1KeEwD8ZamzgYmBHg5wLKjPScQkfOAY4GwRWQu8BfxGRF6vdE4W0AFARBKAxsDOUAbpULVjUdUdqlro/3ICcGRoQ3ROVTf5/9yGry45sNIp+58XvzRgU2iic666cajqblXd6//8YyBRRFqEPNDqZQFZqvqj/+tJ+BJn5XPC/TmpdhwR9JyUOx34RVW3BnjM9eckYhK+qt6pqmmq2gnfr0RfqOoVlU57H7jS//kF/nPCbdbiaCyVandn43txN+yISKqINCz/HDgFWFjptPeB3/u7EAYDuaq6OcShHpKTcYhIm/LXhERkIL7/PztCHWt1VHULsEFEuvsPnQgsrnRa2D8nTsYRKc9JBZcSuJwDIXhOIrFL5wAiMgaYparv43tx5z8ishLfzP4ST4OroUpjuVlEzgZK8I1lhJexHUJr4F3//7kE4L+q+qmI3ACgqv8GPgbOAFYCecBVHsV6KE7GcQHwRxEpAfKBS8JxQuH3f8Ab/hLCauCqCHxOoPpxRMxzIiL1gZOBP1Q4FtLnxJZWMMaYGBExJR1jjDF1YwnfGGNihCV8Y4yJEZbwjTEmRljCN8aYGGEJ30Q0EdnrwjX7i8gZFb6+V0RGVXFuioh87V8Qqy73TBKRb/xvGDTGFZbwjTlYf3z90E5cDUxR1dK63FBVi4AZhPfqribCWcI3UUNEbhWRn/0Lzt3nP9ZJfOuoTxCRRSIyXURS/I8d5T93poiMFZGF/jf4jAEuFt+65eUJuJeIfCUiq0Xk5gq3vRx4r0IMt4lvTf15IvIP/7GvRGScfwa/xH/fKSKyQkQeqHCtqf7rGeMKS/gmKojIKfjWER+Ib4Z+ZIXFz7oBz6hqbyAHON9//GXgBlUdApTC/pn2aHx7EfRX1bf95/YATvVf/x4RSfT/cOiiqmv9MZwOnAsM8i/m9WiFEItUdSjwb3w/IG4E+gAjRKR8RdeFwFHB+jsxpjJL+CZanOL/mINvBcUe+BI9wBpVnev/fDbQSUSaAA1V9Xv/8f9Wc/2PVLXQvznKNnxLMbTA9wOk3EnAy6qaB6CqFRfue9//5wJgkapu9i+Otxr/gln+slBR+Zo+xgSbvUBkooUAD6vqcwcc9G1zWVjhUCmQQuClaA+l8jUS8O23kFwphqrWKin//rJK1yrjwP+H9YCCGsZmjCM2wzfRYhpwtYg0ABCR9iLSqqqTVXUXsMe/KiEcuNDeHqDaWbb/GvEiUp70p/tjqO+PoVlNBuAv7WSranFNvs8Ypyzhm6igqtPxlWVmisgCfGunV5e0rwGeF5GZ+Gbnuf7jX+J7kbbii7ZVmQ4c64/hU3ylm1kiMhcI2Mp5CMPwrZhojCtstUwTs0SkQfnmGSJyB9BWVf9cw2sMAP6qqr8LQjxTgDv9G14bE3RWwzex7EwRuRPf/4N11GLPAVWdIyJfikh8XXrx/R0/Uy3ZGzfZDN8YY2KE1fCNMSZGWMI3xpgYYQnfGGNihCV8Y4yJEZbwjTEmRvw/qeN1vhzB/AUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fun_kernel='linear'\n",
    "\n",
    "X,y = gen_data(data_size=data_size,visualization=visualization)\n",
    "X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2)\n",
    "svm = SVM(max_iter=max_iter, C=C, fun_kernel=fun_kernel, epsilon=epsilon)\n",
    "svm.fit(X_train,y_train)\n",
    "acc = svm.score(X_test,y_test)\n",
    "w1,w2 = svm.weights()\n",
    "b=svm.b\n",
    "print(\"支持向量数为\", len(svm.support_vectors)) \n",
    "print(\"svm预测准确率：\",acc)\n",
    "print(\"w权重分别表示为\", w1,w2,b)\n",
    "print(\"可视化分类界面如下\")\n",
    "show_result(X,w1=w1,w2=w2,b=b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "svm预测准确率： 1.0\n"
     ]
    }
   ],
   "source": [
    "# 测试多项式核的结果\n",
    "fun_kernel='linear'\n",
    "svm = SVM(max_iter=max_iter, C=C, fun_kernel=fun_kernel, epsilon=epsilon)\n",
    "svm.fit(X_train,y_train)\n",
    "acc = svm.score(X_test,y_test)\n",
    "w1,w2 = svm.weights()\n",
    "b=svm.b \n",
    "print(\"svm预测准确率：\",acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "svm预测准确率： 0.9\n"
     ]
    }
   ],
   "source": [
    "# 测试多项式核的结果\n",
    "fun_kernel='gaussian'\n",
    "svm = SVM(max_iter=max_iter, C=C, fun_kernel=fun_kernel, epsilon=epsilon)\n",
    "svm.fit(X_train,y_train)\n",
    "acc = svm.score(X_test,y_test)\n",
    "w1,w2 = svm.weights()\n",
    "b=svm.b \n",
    "print(\"svm预测准确率：\",acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "svm预测准确率： 0.954\n"
     ]
    }
   ],
   "source": [
    "# 测试高斯核在mnist数据集上的表现\n",
    "# 这里使用mnist的前1000个数据作为训练数据\n",
    "max_iter = 150\n",
    "fun_kernel='gaussian'\n",
    "X_train,X_test,y_train,y_test = load_mnist(\"./data/mnist_train.csv\",\"./data/mnist_test.csv\")\n",
    "svm = SVM(max_iter=max_iter, C=C, fun_kernel=fun_kernel, epsilon=epsilon)\n",
    "svm.fit(X_train[:1000],y_train[:1000])\n",
    "acc = svm.score(X_test[:1000],y_test[:1000])\n",
    "print(\"svm预测准确率：\",acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
